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ABSTRACT 

Radiative cooling by metals in shocked gas mediates the formation of ionization front instabilities 
in the galaxy today that are responsible for a variety of phenomena in the interstellar medium, 
from the morphologies of nebulae to triggered star formation in molecular clouds. An important 
question in early reionization and chemical enrichment of the intergalactic medium is whether such 
instabilities arose in the H II regions of the first stars and primeval galaxies, which were devoid of 
metals. We present three-dimensional numerical simulations that reveal both shadow and thin-shell 
instabilities readily formed in primordial gas. We find that the hard UV spectra of Population III 
stars broadened primordial ionization fronts, causing II2 formation capable of inciting violent thin- 
shell instabilities in D-type fronts, even in the presence of intense Lyman- Werner flux. The high post- 
front gas temperatures associated with He ionization sustained and exacerbated shadow instabilities, 
unaided by molecular hydrogen cooling. Our models indicate that metals eclipsed II2 cooling in I-front 
instabilities at modest concentrations, from 1 x 10"'^ - 1 x 10~^ 7jq. We conclude that ionization 
front instabilities were prominent in the H II regions of the first stars and galaxies, influencing the 
escape of ionizing radiation and metals into the early universe. 

Subject headings: H II regions: simulation — cosmology: theory — early universe 



1. INTRODUCTION 

Instabilities of several varieties are known to form 
in ionization fronts (I-fronts) in the galaxy today, 
giving rise to structures like elephant trunks and 
bright rims in star-forming regions of the interstellar 
medium (ISM) . Two prom inent types are shadow in- 
stabilities (Williams 1999) and thin-shell instabilities 
(jCarcia-Segura & Franco 1996), which usually manifest 
as fingers of ionized gas protruding ahead of the front. 
Shadow instabilities form when a density fluctuation is 
advected through an R-type front, creating a dimple 
that first elongates but then stabilizes. When the front 
becomes D-type, the dimple can erupt in a jet if the 
shocked gas radiativcly cools and forms a cold dense shell 
liable to fragmentation. Radiation escapes through fis- 
sures in the shell, causing the dimple to fracture and 
jet forward. Thin-shell instabilities arise when a D-type 
front encounters a density perturbation. If the prefront 
shocked gas efhcientl y cools, t he perturbation triggers a 
Vishniac instability ( Vishnia3 [l983,) in the shell; if not, 
the shell may instead fracture due to R ayleigh- Taylor in- 
stabilities ()Mac Low fc Normanlll993h . Either way, radi- 
ation rapidly destabilizes the perturbations by escaping 
through cracks in the shell and driving ionized exten- 
sions of gas much larger than the original disturbance 
into the surrounding medium. In a previous paper, we 
examine d the topology of both ty pes of instabilities in 
the ISM (jWhalen fc Normanll2007n . confirming that their 
violence was determined by the cooling efficiency of the 
shocked gas. 

H II regions a ppeared at the redsh i fts of first struc- 
ture formation (jWhalen et al. I l2004t iKitavama et al. I 
[20nl . be ginning with the formation of isolated mas- 
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sive stars in the first dark ma tter halos of mass > 1.0 
X IQS M0 at z 7 20 - 30 (lAheLeTaD [2OOOI I200I : 
iBromm et al. I[l999l : iBromm et al.ll2002( ) and continuing 
with the rise of the ffi'st stellar populations in primeval 
galaxies a few hundred Myr later. Whether I-front in- 
stabilities could form in pristine H and He at high red- 
shifts is a key question in early cosmological reioniza- 
tion and structure formation for several reasons. First, 
they may have regulated UV escape from the first stars 
and galaxies, which has risen sharply in importance with 
the recent discovery of substantial free electron opti- 
cal depths at by the Wilkinson Microwave Anisotropy 
Probe (WM AP) (lSDergel"ei^[200a iKogut et al.ll2008l : 
iPage et aD l2006l ). implying that cosmological reioniza- 
tion began much earlier than previously expected. Sec- 
ond, I-front instabilities may have created density con- 
densations prone to clumping and later collapse, forming 
new generations of stars. Third, the mixing of metals 
from the first supernovae with remnants of I-front in- 
stabilities in relic H II regions determined the extent to 
which the early intergalactic medium (IGM) was chem- 
ically enriched. Contamination of instability relics may 
fracture and collapse them on ma ss scales very differen t 
from the stars that created them iMackev et al.l (|2003f l: 
lYoshida et al.l ()2007f ). abruptly splitting the primordial 
initial mass function (IMF) into a low-mass branch in 
just one generation. On the other hand, if saturated I- 
front instabilities drive turbulent flows, they may have 
supported clouds in minihalos or protogalaxies from col- 
lapse , postpon ing rather than promoting star forma- 
tion (IWhal en fc Norman 2007). Finally, recent work by 
iMizuta etal, (2005 ., _20_Q6) suggests that ionization fronts 
penetrating molecular cloud cores in the ISM may de- 
velop unstable modes that reach the centers of the core 
before the front. This may disrupt collapse of cloud cores 
into new stars, both in the galaxy and at high redshifts. 

Given the central role of radiative cooling in modulat- 
ing I-front instabilities in the ISM, one might question 
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whether they arose at aU in primordial ionized flows, in 
which the atomic or molecular hydrogen cooling is far 
less efficient. Cosmological ionization fronts also differed 
from their galactic counterparts in being driven by the 
hard UV and soft X-ray spectra of Pop 111 stars and 
miniquasars. Such spectra exhibit a range of photon 
mean free paths that can, depending on the ambient den- 
sity, greatly broaden the front, in some cases to as much 
as the Hubble distance. It is unclear what impact the 
broadening of the front would have on the formation and 
evolution of unstable modes in I-fronts; the temperature 
precursor preceding such fronts may preempt their for- 
mation entirely. 

We have performed three-dimensional radiation hy- 
drodynamical calculations of ionization front instabili- 
ties in which multifrequency radiative transfer is cou- 
pled to nine species primordial chemistry. We explore 
the formation of thin-shell and shadow instabilities in 
both primordial and incrementally-enriched gas in ide- 
alized geometries to disentangle their evolution from 
other features imprinted on the front by density imho- 
mogeneities in real cosmological objects. As primor- 
dial fronts exit minihalos and protogalaxies they will en- 
counter structures on several spatial scales that shape 
the global morphology of the final H II region (e.g., 
lAbel et all ([2006)). For example, halos hosting Popu- 
lation III (Pop III) star formation exhibit filamentary 
inflows an d other dep artures from spheric ity on 100 
PC scales (lAbel et al. \ [2 000; Fuller & Couc hmanI 120001 : 
lYoshida. Abel. Hernquist. fc Sugiyamaii2003i) . However, 
protostellar infall likely settles into an accretion disk 
on 1000 AU scales, collimating the emergent front into 
biconical flows above and below the plane of the disk 
(|Tan &: McKed 120041 ) . Transverse gas velocities inher- 
ited from the angular momentum of the parent dark mat- 
ter halo further complicate breakout of the nascent front 
on intermediate scales. We survey the formation of in- 
stabilities in primordial fronts in advance of multiscale 
simulations of UV breakout from primeval stars and and 
galaxies. 

In § 2 we describe the multifrequency upgrades 
to our photon-conserving radiative transfer in the 
ZEUS-MP reactive flow physics code documented in 
IWhalen fc Normal] (|2006l ). Primordial cooling processes 
are examined in § 3, particularly the role of H2 cool- 
ing in cosmological ionization fronts. In § 4 we present 
our numerical models of D-type I-front instabilities due 
to thin-shell modes with atomic and molecular cooling. 
We follow the evolution of shadow instabilities associ- 
ated with Pop III blackbody spectra in § 5, pinpoint the 
threshold at which metal line cooling triggers thin-shell 
instabilities in D-type fronts in § 6, and conclude in § 7. 

2. RADIATIVE TRANSFER: CONSERVATIVE 
MULTIFREQUENCY PHOTON TRANSPORT 

The monoenergetic radiative transfer in our previous 
studies does not capture postfront gas temperatures from 
first principles or the natural broadening of fronts by 
hard UV photons. It also neglects molecular hydro- 
gen dissociation by Lyman- Werner photons and photode- 
tachment of H^, key processes that regulate the abun- 
dance of H2, an important coolant in primordial gas. 
Here we describe modifica tions to our numerical method 
(jWhalen fc NormanI |2006[) for computing radiative rate 



coefficients with a multifrequency algorithm that pre- 
serves photon conservation. 

Constraining the number of absorptions in a zone due 
to all processes to equal the number of photons entering 
the cell minus those exiting leads to 

Uabs oc 1 - e~'^, (1) 

where 

n 

T^J2 '^^^^'^r, (2) 

i=l 

is the total optical depth of the cell (the sum is over 
all removal processes). Equations 1 and 2 are true for 
a single energy and therefore hold if summed over all 
energies. The total absorption rates riabs are necessary 
to determine the flux that reaches the next zone. On 
the other hand, the rate for a single removal process is 
a function of an attenuation factor that does not involve 
any of the other types of absorptions: 

rii (X 1 — e^"^*, (3) 

where = nidiAr, the optical depth to that interaction. 

The sum of the individual rates evaluated in this man- 
ner is inconsistent with Hats, but the two can be recon- 
ciled by enforcing 

ni tX = Uabs- (4) 

J2^=l 1 - e 

This formulation guarantees that the sum of the rates 
matches the total rate associated with the optical depth 
of the zone (thus preserving photon conservation) while 
properly binning individual reactions according to their 
numbers. We evaluate rate coefficients (where k^ = 
Hi / (nV cell) , where n is the gas number density) at ev- 
ery frequency, assuming the photon emission rate n{i>i) 
at that frequency, and sum them to obtain the total 
coefficient for the entire spectrum. We increment en- 
ergy release in the gas one frequency at a time for a 
given photoreaction as well. In tests spanning 40 to 2000 
logarithmically-spaced bins with a 1 X 10^ K Pop III 
blackbody spectrum, we obtain good convergence with 
80 bins above the ionization limit of H. 

We partition the emissi on rate h-y of ion izing photons of 
a Pop III star above from lSchaereil (|2002f ) into frequency 
intervals according to the number rate of a blackbody: 

4=1 i=l ^ ^ 

where the sum is over \w > 13.6 eV. The emission rate 
for a particular bin is then 

n^j = ri(z^j)Ai/ = n^. (6) 

The procedure is similar for a quasar power-law flux spec- 
trum F(z^)dz^ oc v^^dv. If fiQso is the ionizing emission 
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TABLE 1 

ZEUS-MP Radiative Reaction Processes 



Energy 



Rate 




Reaction 


k24 


H + 


7 


H+ + e- 


k26 


He+ 


+ 7 


He++ + e 


k2S 


He - 


f 7 - 


* He+ + e- 


k27 


H- 


+ 7 - 


H + e- 


k28 




+ 7 - 


~» H+ + H 


k29 


Ha - 


f 7 - 


H+ + e- 


kso 




+ 7 - 


2H+ + e- 


kai 


Ha - 


f 7 - 


-* 2H 



hi/ > 13.6 eV 

hu > 54.4 eV 

hu > 24.6 eV 

hu > 0.755 eV 

2.65 eV < h!^ < 21 eV 

hu > 15.42 eV 

30 eV <hu < 70 eV 

11.18 eV <hu < 13.6 eV 



rate for the quasar, a simple integration over frequency 
yields 



from which follows 



Q + l 



nqso 



,,.a + l 



(7) 



(8) 



We truncate the ionizing UV spectrum at 90 eV for Pop 
III stars. 

The radiative processes now in the code are listed in 
Table [H We adopted cross sections from lOsterbrockl 
(|1989f) to calculate k24 and k26, from 'Abel et all ()1997f ) 



to compute k27 and from Shapiro & Kang ( 1987|) to eval- 
uate k25, k28, k29, and kao. The rate coefficients k27 and 
k28 both extend below 13.6 eV; we allocate an additional 
40 bins equally spaced between 0.755 eV and 13.6 eV to 
evaluate them. We obtain binned source rates over this 
energy range by constructing blackbody number spec- 
tra, normalizing them by total stellar luminosities taken 
from ISchaereij ( 20021 ) instead of by total ionizing pho- 
ton emission rates as before. Considerable speedups of 
the radiative transfer 30%) can be realized if the lin- 
ear approximations to equations [1] and [3] are adopted 
in the limit of small optical depths, which arc common 
for species like and HJ. The kai calculation is de- 
coupled from the radiative transfer and computed us- 
ing th e H2 self-shielding functions of iDraine fc Bertoldil 
(Il996t) corrected for thermal Doppler broadening (see 
also I Abel et"all (|1997[ )) 



fcgi = 1.1 X 10^F{l^)F, 



shield •) 



(9) 



where F shield is obtained from eq. 37 of lDraine fc Bertoldil 
(|1996[ ) and F{i^) is the specific fiux in s~^ ergs cm~^ Hz~^ 



(10) 



and the frequency range corresponds to the Lyman- 
Werner band 11.18 eV to 13.6 eV. Total stellar luminosi- 
ties are applied to derive F(z^) for the blackbody spec- 
trum. To approximate the effect of gas motion we assign 
a temperature of 10,000 K to the Doppler correct ion to 
the shielding function, as in lAhn fc Shapiro! (l2007h . The 
viability of thermal broadening as a proxy for flows in 
Tfronts is uncertain, but including it reduces shielding 
at intermediate column densities, setting lower limits on 
H2 cooling in the gas. 
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Fig. 1. — Primordial H and He abundances in a cosmological 
ionization front at t = 0.6t« . In this figure, the radiation wave 
approaches from the right towar d the halo centered at the left 
boundary. Compare to Fig. 8 of 'Ahn S hapirol II2007I ) and the 
mirror imag e of Fig. 3 of lRicotti"et al. 1,2001,) . Each curve is labeled 
by its species. 



To test the coupling of the multi frequency tra nsport to 
the nine-species reaction network (jAnninos et al. 1997), 
we simulated the photoevaporation of a cosmological 
minihalo 540 pc from a 120 M© Pop III star . A 
spherically-symmetric truncated isothermal sphere (TIS) 
density profile was centered at the left x face of a 125 
pc cartesian box with a uniform mesh of 512 zones. A 
geometrically attenuated plane wave propagated along 
the negative x-axis, engulfi ng the halo. Th i s is th e test 
case shown in Fig. 8 of lAhn fc Shapirol (|2007f ). and 
we adopted identical initial conditions for comparison 
to their one-dimensional lagrangian radiation hydrocode. 
Number fraction abundances for the primordial H and He 
species appear in Fig. [l]at 0.6t*, where t* is the lifetime 
of the star. 

Perhaps the most striking feature is the shell of molec- 
ular hydrog en created in the pa rtially ionized outer layers 
of the front (jRicotti et al.ll200r ) (notice its peak coincides 
with an H II fraction ~ 0.1). The hard UV photons in 
the spectrum transit a range of mean free paths in the 
neutral gas, broadening the front to approximately 5 pc. 
Free electron abundances at temperatures of a few thou- 
sand K in the finite width of the front rapidly catalyze 
II2 production by the and channels: 



H + e~ ^ H~ + J 

H + H+ ^H++j 



H-+H^H2 + e- (11) 
H+ +H^H2 + H+ (12) 



The peak and abundances are in step with the 
II2 peak, as expected since they are the key interme- 
diaries for molecular hydrogen production in these rela- 
tively low densities. The channel is dominant because 
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radius (pc) radius (pc) radius (pc) 

Fig. 2. — Left: Primordial cooling rates in the ionized gas and neutral shocked shell of an I-front at 100 Kyr not cooled by H2. Dalgarno 
& McCray cooling rates, included for comparison, dwarf atomic H radiative rates in the dense shell. Right: H2 cooling rate, density, 
temperature, and velocity profiles for a primordial ionization front with no Lyman- Werner photodissociation (solid) and Lyman- Werner 
photodissociation (dotted) . 



it is more abundant and it has much higher H2 forma- 
tion rates in th ese regimes. Fig . [T] is in good agreement 
with Fig. 8 of I Ahn fc Shapiro! (|2007D and is similar to 
Fig. 3 of iRicotti et all (|200H ). which depicts H and He 
abundances for a primordial ionization front in a static 
uniform medium at the mean cosmological density (their 
figure should be compared to the mirror image of ours 
because their front propagates to the right). We note the 
relative abundances of and to be a general fea- 
ture of partially ionized primordial gas at a few thousand 
K at t hese densities (compare Fig. 3 in Anninos et oFI 
(|1997| ) , which features coUisionally ionized gas collapsing 
in a pancake shock). 

3. SHOCK COOLING IN PRIMORDIAL GAS 

Five mechanisms that cool primordial H and He in 
I-front shocks are electron coUisional and ionizational 
cooling, recombinational cooling, bremsstrahlung, and 
Compton cooling (primarily in the early universe). The 
free electrons required for these processes are present in 
minute quantities (ue/u// ~ 1 x 10"'') at high redshifts 
and are either remnants of the recombination era or cre- 
ated in cosmological shocks. These channels cool H and 
He far less efficiently than metal lines in the ISM at T 
< 5000 K. However, as discussed in § [2l molecular hy- 
drogen may f orm in the I-fronts of very massive Pop HI 
stars (jRicotti et al.ll2001[ ) and significantly raise cooling 
rates in the neutral gas. H2 ro-vibrational lines cannot 
rival metal cooling (since they can only lower primordial 
gas to ^ 200 K instead of 10 - 20 K for metals) but are 
much more effective than H and He alone at temperatures 
below 5000 K. H2 formation is greatest in the partially 
ionized gas at the base of the shocked shell where energy 
losses may collapse the shell and trigger dynamical per- 
turbations. However, the role of H2 cooling in primordial 
I-front instabilities is unclear because massive stars with 
strong ionizing luminosities are also intense sources of 
Lyman- Werner (LW) photons that dissociate H2. With 
energies below the ionization limit, these photons can 



cross the front and degrade the molecular hydrogen cat- 
alyzed in its outer layers. 

We formulated three sets of one-dimensional calcula- 
tions to examine I-front shocks when H, He and H2 cool- 
ing are present and absent. Collapse of the shell in one 
dimension is a rough predictor of instability formation in 
three dimensions. The initial density profile was: 



n(r) 



Tic if r < Tc 

nc(r jr^)^'^ if r > 



with Uc and equal to 1 x 10"* cm~^ and 0.2 pc, re- 
spectively. 

A UV source with h^h = 1-0 x lO"'* s^^ was centered in 
the mesh, ensuring the front would become D-type within 
the core. The grid was discretized into 200 uniform radial 
zones with inner and outer boundaries of 0.01 pc and 
2.0 pc and reflecting and outflow boundary conditions, 
respectively. The gas was set to 500 K, appropriate to 
star-for ming mi nihalos cooled by molecular hydrogen at 
z ~ 20 ()Abel e t al. 2000, 2002). We chose a cooling 
cutoff (defined to be the gas temperature below which 
cooling curves were not applied to gas energy updates) 
equal to the cosmic microwave background Tcmb =2.73 
(l+z) K with z = 20. 

All three runs were multifrequency, with black- 
body spectra and nin e-species primordial gas chemistry 
([Anninos et "anil997[ ). The gas profile was initialized 
with H and He number densities whose sum was equal to 
Uc but 76% and 24% H and He by mass, respectively. We 
adopted initial electron and H2 number fr actions of 1 x 
10"^ and 2 x 10"'^ consistent with 2-20 ()Ricotti et al.l 
200I. When H2 cooling was included, the curves of 
Galli fc Fallal (fl998h were utilized. The ionizing photon 
emission rate corresponds approximately to a star with a 
surface temperature, mass, and lum inosity of 57,000 K, 
15 M© and 2.1 x lO'' L© (|Schaereil |20l)3l respectively 
Wc applied these parameters to normalize the photon 
rates binned in each frequency interval. 



Fig. 3. — Evolution of a D-type front in an primordial r ^ envelope. Panels (a) and (b): temperatures at 47.6 kyr and 135 kyr, 
respectively. Panels (c) and (d): densities at 47.6 kyr and 135 kyr, respectively. 



3.1. H/He Cooling in Ionized and Shocked Gas 

H2 was allowed to form but not cool the gas in the 
first run. As an experiment, we toggled the five atomic 
H and He cooling processes to determine their effect on 
the temperature of the ionized gas. We found (as ex- 
pected) recombinational cooling to be dominant, reduc- 
ing the postfront temperature from 60,000 K to 30,000 
K. The remaining four channels collectively cooled the 
gas another 4000 K, down to 26,000 K. AU five rates 
are shown in in Fig. [2] at 100 Kyr in the ionized and 
shocked neutral gas. The 1-front ends at 0.94 pc in the 
plot and the shocked shell extends to f .19 pc. While coUi- 
sional ionizations briefly surpass recombinations near the 
front, bremsstrahlung, interestingly, is the second most 
important cooling mode in most the of the H II region, 
primarily due to the relatively large temperatures associ- 
ated with He ionization. Collisional excitation of neutral 
hydrogen then follows collisional ionization in order of 
importance in the postfront gas, with Compton cooling 
as the least effective. 

Bremsstrahlung and collisional ionization and exci- 
tation drop sharply in the shocked gas, leaving only 
Compton and recombinational cooling to operate in 
the low free electron fractions ther e. For the s ake o f 
comparison, we overlay the IDalgarno fc McCravl ()1972l ) 



cooling rates for an elec tron fraction of 0.01 used in 
IWhalen fc NormanI ()2007[ ). Metal line cooling is eight or- 
ders of magnitude greater than recombinations or Comp- 
ton scattering. Indeed, if no cooling cutoff whatsoever 
was employed with the metal lines they would quickly 
collapse the shocked gas to very dense cold shell whose 
final thickness is bounded from below only by numeri- 
cal resolution. We note that the Compton component is 
highly dependent on redshift; had this simulation been 
performed at z = 35, inverse Compton scattering would 
have dominated recombinations in the postshock shell. 

Recombinational and Compton cooling are so low in 
the shocked gas that they are inconsequential to its 
structure. The density and temperature profiles of the 
shell are indistinguishable from those of a nonradiating 
shock. Instability formation, except possibly those due 
to Rayleigh- Taylor perturbations, is extremely unlikely 
in these circumstances. 

3.2. H2 Cooling 

In the second model H, He and molecular hydrogen 
cooling is included but not Lyman- Werner dissociation, 
so H2 forms freely in the outer layers of the front. In 
the third model we activate H2 photdissociation as de- 
scribed in § [2I H2 cooling rates, temperatures, densities. 
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Fig. 4. — Evolution of a D-type front in an primordial r envelope. Panels (a) and (b): temperatures at 47.6 kyr and 123.7 kyr, 
respectively. Panels (c) and (d); densities at 47.6 kyr and 123.7 kyr, respectively. 



and velocities for both runs appear in Fig. [21 Molecular 
hydrogen efficiently drops the base of the shocked shell to 
700 K in the absence of dissociating radiation, condens- 
ing the shell to higher peak densities than in the first set. 
It is likely in these circumstances that H2 cooling would 
cause unstable modes and shell breakup in primordial 
I- fronts. 

However, at this proximity to the central star molec- 
ular hydrogen does not fully self-shield from the dissoci- 
ating UV radiation, as shown in our third model. The 
sharp reduction in cooling rates attests to the efficiency 
with which LW photons destroy the H2 molecules (low- 
ering their abundance by more than two orders of mag- 
nitude), resulting in a shell structure indistinguishable 
from one with atomic H cooling only. Comparison of 
the density profiles indicates that H2 cooling reduces the 
shell thickness by half in the absence of LW radiation and 
that it approaches the efficiency the Dalgarno & McCray 
curves within a thin layer at the base. The velocity pro- 
files and positions of the two fronts indicates that H2 
cooling slows the front by 10% in the absence of LW 
flux. Nevertheless, even with photodissociation the H2 
cooling rates are comparable to the Dalgarno & McCray 
rates at electron fractions of 0.0001 - 0.001, suggesting 
that even low molecular hydrogen abundances set by LW 
fluxes may induce unstable modes in primordial Lfronts 



(jWhalen fc NormanI [200I . an issue to which we return 
below. 

4. DYNAMICAL INSTABILITIES IN D-TYPE PRIMORDIAL 
I-FRONTS 

Can atomic cooling in H and He alone incite dynam- 
ical instabilities in D-type primordial ionization fronts? 
In this section we address both Vishniac and Rayleigh- 
Taylor instabilities driven by ionization fronts in per- 
turbed spherically-symmetric density gradients. 

4.1. Vishniac Modes 

Our one-dimensional model in § 13.11 revealed that I- 
front shocks cooled only by atomic H and He were indis- 
tinguishable from nonradiating shocks and therefore are 
not susceptible to Vishniac instabilities. To verify this 
prediction we performed a three-dimensional simulation 
with monoenergetic photons to prevent any complica- 
tions due to multifrequency broadening of the front. In 
this run the same central photon rate, radial mesh and 
density profile were used as before but with 180 zones in 
theta and in phi, with reflecting boundaries at 7r/4 and 
37r/4 in both angles. The initial gas and cooling cutoff 
temperatures were set to the CMB background at z = 20 
and to 1000 K, respectively. The energy per ionization 
er was chosen to be 1.6 eV in order to set the postfront 
temperature to 10000 K. Thus, while He was present in 



Fig. 5. — Temperature images of H2 mediated thin-shell instabilities in primordial ionization fronts. Left: no LW radiation, 40 kyr. 
Right: LW radiation included, 51.4 kyr. 
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Fig. 6. — IGalli &: Fallal I I1998I ) molecular hydrogen cooling rates 
for three fiducial densities. Solid: n//2 = 1 x 10^ cm~^; dotted: 
= 1 X 10^ cm~^; dashed: ng^ = 1 x 10** cm""^. 



the calculations to preserve its inertia to hydrodynamical 
motion, it was not ionized. Ignoring He ionization sim- 
ply results in lower ionized gas temperatures, which do 
not deter instabilities from forming (He contributes even 
less to I- front shock cooling than H in these regimes). 

To seed the formation of instabilities, we randomly var- 
ied the gas density of each cell by at most 1%, holding gas 



energies constant to prevent pressure fluctuations from 
prematurely smoothing the perturbations. We only im- 
posed these variations beyond radii of 0.125 pc to prevent 
the onset of shadowing instabilities in the R-type front 
before dynamical instabilities in the D-type front could 
form. In Fig. [Sj the temperatures and densities reveal 
small fluctuations in the ionization and shock fronts due 
to the original perturbations, but they do not grow with 
time. The shell simply becomes thicker as more shocked 
neutral gas accumulates. We note the layer of gas pho- 
toevaporating into the interior of the H II region from 
the inner face of the shell in panel (a) of Fig. \3\ It is 
visible as the thin red arc of lower densities bounding 
the yellow layer of denser shocked gas from below and 
can only be seen in the absence of instabilities. Our cal- 
culation confirms that dynamical instabilities in I-fronts 
cannot arise through the Vishniac instability in shocked 
gas cooled by only by H and He. Since our result de- 
pends only on cooling efficiencies, they will hold for the 
wide range of densities encountered in primordial halos 
and protogalaxies. 

4.2. Rayleigh- Taylor Modes 

Rayleigh- Taylor (RT) instabilities can grow in primor- 
dial I-fronts if they accelerate in a density gradient, which 
is the case for profiles steeper than r~^. We performed a 
run with parameters identical to those in i) l4.1l except in a 
density gradient of r~^'^ and show (f) = tt/A temperature 
and density slices of its evolution in Fig. [H Perturba- 
tions in the dense shell are somewhat more pronounced 
than in the Vishniac run, hinting at the onset of unstable 
fingers, but they never fully materialize and disrupt the 
shock as the front exits the envelope. The shock clearly 
gains strength down the gradient, reaching temperatures 
of nearly 10,000 K. In such regimes electron-neutral col- 
lisions soon ionize the shock, helping the front to break 
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Fig. 7. — Density images of a shadow instability due to an underdense perturbation in atomic H with hv = 14.2 eV and no H2 cooUng: 
(a) 1.9 kyr, (b) 4.9 kyr, and (c) 12.6 kyr. 



through. While possible in general, we find that RT in- 
stabilities fail to induce unstable modes in the front at 
these densities, which are typical of those in 1 x 10^ - 
1 X 10^ M0 minihalos at early epochs and of molecular 
cloud cores in protogalaxies at later times. 

4.3. Thin-Shell Modes Due to H2 Cooling 

We show in Fig. Othe evolution of thin-shell instabili- 
ties in primordial D-type ionization fronts due to molec- 
ular hydrogen catalyzed at the base of the shocked shell. 
The first run was identical to that of § 14.11 except that 
H2 cooling was activated, but not LW photodissociation. 
The second was identical except Lyman- Werner flux was 
included. These models span the range of self-shielding 
anticipated in cosmological halos and protogalaxies, from 
virtually none to complete. In the absence of LW pho- 



tons, H2 abundances of 1 x 10"'^ to 1 x 10^^ form at the 
Lfront/shock interface and, as shown in the left panel of 
Fig. [5l are as effective as metals at instigating dynam- 
ical instab ilities in the D-type front (compare with Fig. 
5 of Whalen &: Normanl ([20071 )). When LW flux is in- 
cluded, equilibrium H2 fractions in the shell fall to 1 x 
10~^ to 1 X 10~^ but, surprisingly, strong instabilities 
persist despite the diminished cooling. The key to this 
phenomenon is the temperature dependence of the H2 
cooling curves, shown in Fig. [6l From 300 K, the typical 
temperature of a star-forming cosmological minihalo, to 
3000 K, the temperature of the shock in the primordial 
front, H2 cooling rates rise by three orders of magnitude, 
offsetting losses due to the photodissociation of the shell. 
We find that the instabilities are not as violent but still 
heavily perturb the shock. 



Fig. 8. — Shadow instability due to an underdense perturbation, full nine species primordial chemistry with a 100,000 K blackbody 
spectrum: (a) 1.3 kyr, (b) 6.3 kyr, and (c) 13.4 kyr. 



Since our assumed central densities are lower than 
those in which Pop III stars form, we expect that H2 
will be more heavily shielded and thin-shell instabilities 
will be more violent in primordial halos and protogalax- 
ies than in this second case. Preliminary calculations in 
1 X 10^ Mq halo profiles hosting 120 stars confirm 
this suspicion, even with their much larger LW fluxes. 
These two cases likely bracket the magnitudes of the in- 
stabilities that occurred in primordial D-type ionization 
fronts. 

5. SHADOW INSTABILITIES IN R-TYPE PRIMORDIAL 
I-FRONTS 

To test if the shadow instability can form in primordial, 
we performed a calculation in which a plane-parallel R- 
type ionization front enters the yz-face of a rectangular 



box with a uniform density except for a spherical un- 
dcrdcnsity that is slightly offset from the face of entry. 
Its density decreases linearly in radius from the ambient 
value at its center to 50% below this amount at its sur- 
face. Our simulation mesh is 1.0 pc along the x-axis and 
0.25 pc along the y and z-axes, with 1000, 250, and 250 
zones in the x, y, and z directions respectively. The gas 
density in the box is 1000 cm^^ with a temperature of 72 
K to establish a sound speed c^ of 1 km/s. The incident 
photon flux along the x-axis at the face of entry is 3.0 
X 10^^ cm~^ s~^ to approximate that of a 120 Mq star 
at 1.71 pc. The radius of the perturbation is 0.0125 pc 
and is placed 0.05 pc along the x-axis and centered in the 
yz-plane. Inflow and outflow boundary conditions were 
assigned to the pc and 1.0 pc faces, respectively, with 
reflecting boundaries on the other four faces. 
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5.1. Monoenergetic Spectrum/ Atomic Cooling 

We first consider shadow instabilities in H gas with 
monoenergetic photons to avoid front broadening, with 
fixed energy per ionization ep = 0.8 eV to ensure post- 
front temperatures of 10000 K. We show xy density shces 
at z = of the instability in Fig. [7]at 1.9 kyr, 4.9 kyr, and 
12.6 kyr. By 1.9 kyr the front has become D-type and 
the original corrugation in the R-type front is now non- 
linear. The expanding photoevaporated perturbation is 
faintly visible on the left in the figure. The double peaks 
are cross sections of a ringed jet whose morphology is due 
to the radial density profile of the sphere. The R-type 
front preferentially advances along lines of sight parallel 
to the X-axis that cut the underdense regions close to the 
surface of the sphere. Along the axis piercing the sphere 
through its center, the front advances at nearly the same 
rate as in the unperturbed medium because the densities 
along this path through the sphere are closer to those in 
the surrounding gas. 

The evolution of the insta bility is qualitatively simila r 
at first to that in Fig. 8 of IWhalen fc NormanI (|2007f ). 
but it quickly dampens and becomes subsumed in the in- 
efficiently cooled shocked shell. The density knots in the 
shock located above and below the midp l ane o f the jets 
at their base in the IWhalen fc Normal] ()2007[ ) run are 
also visible in this run at first but they rapidly dissipate 
with little lateral migration a cross the face of the shock. 
If the assertion by 'William^ (|1999( ) that these knots are 
odd-even numerical instabilities is true, they cannot sur- 
vive long in the steadily thickening shocked layer at these 
densities. Rather than condensing and fragmenting, the 
shell tends to expand back into the jet, crowding its rar- 
efied interior with gas that shields its tip from incident 
radiation and collapsing it back into the face of the shock. 
The collapse of the jet back into a dimple imparts con- 
siderable vorticity to the shocked gas in the vicinity of 
the disturbance, as is evident in panel (c). We find that 
atomic hydrogen cooling effectively quenches the shadow 
instability, but not before it significantly churns the face 
of the shock. 

5.2. Blackbody Spectrum/ H2 Cooling 

The jet instability is far more violent when full nine- 
species primordial chemistry, the blackbody spectrum of 
the 100,000 K primordial star, and molecular hydrogen 
cooling are included, as we show in Fig. [8l Here, it is 
the higher post-front gas temperatures associated with 
multifrequency ionization of H and He that sustain and 
destabilize the jet. At 35,000 K, the high pressure of the 
ionized gas prevents the shocked shell 10,000 K) from 
expanding and backfilling the rarefied jet, so the dimple 
never heals or collapses back into the face of the shock. 
The role of H2 cooling is minimal. The thickness of the 
shell in this model is indistinguishable from that with 
atomic cooling only, but the persistence of some molecu- 
lar hydrogen likely promotes some small-scale fragmen- 
tation of the shell. We observe no saturation of the jet 
and it continues to grow with time. Maximum overdensi- 
ties due to clumping of gas are approximately 10 in both 
runs. 

6. METALLICITY THRESHOLD FOR D-TYPE 
INSTABILITY FORMATION 



In § 4 it was shown that violent dynamical instabili- 
ties do not manifest in shocked shells cooled by atomic H 
only. At what metallicity do dynamical instabilities ap- 
pear in D-type I-fronts in the absence of H2 cooling? We 
performed a grid of calculations with Dalgarno & Mc- 
Cray cooling in metallicities of 1 x 10"** Z0, 1 x 10~^ 
Z0, 1 X 10~^ Z0, 0.1 Z0, and 1.0 Z0 to ascertain if there 
is threshold to instability growth. The problem setup is 
the same as in the three-dimensional r~^ run described 
in § 14. 1[ again assuming an initial electron fraction of 1 
X 10^** appropriate for z ~ 20. At such low fractions 
the cooling in the gas is chiefly due to excitation of met- 
als by collisions with neut ral hydrogen (see Fig. 2 in 
iDalgarno fc McCravl (|1972f )). At solar metallicity these 
cooling rates lie approximately an order of magnitude 
below the Dalgarno fc McCray curve shown in Fig. [2l 
Our metallicities were chosen to survey cooling efficien- 
cies between the very low purely atomic rates and quite 
large 1 Z0 metal line curve. We apply a cooling cutoff 
temperature of 1000 K rather than the CMB tempera- 
ture. This limits the thickness to which the shell can 
collapse so the metallicity at which instabilities appear 
is truly the threshold for their formation. 

In Fig. [5] the I-front and shock are shown at 125.5 
kyr for our six choices of metallicity. We include in 
panel (a) the front in a zero-metallicity g refer- 
ence case in which no unstable modes arise. Noticeable 
perturbations appear in the front in roughly equal num- 
bers in gas enriched to 1 x 10^** Z0 and 1 x 10^'^ Z© 
(panels (b) and (c), respectively); those in panel (c) are 
slightly more pronounced but in neither case do they 
grow significantly with time or penetrate the shell to dis- 
rupt the shock. The agitation in the front does perturb 
the shock in the 1 x 10^^ Z© gas, with visible hydro- 
dynamical motions above and below the image plane. 
Also, consolidation of initially short-wavelength modes 
into fewer larger- wavelength structures is evident, clearly 
de monstrating the evolution in D-type fronts predicted 
bv lGiuHanil (|1979[) . These features are even more pro- 
nounced in the 0.1 Z© case, and the instability com- 
pletely deforms the shock in the 1.0 Z© gas. 

If one adopts deformation of the shock and wavelength 
evolution of the unstable modes as the criteria for in- 
stability formation, it appears that metallicities of ~ 1 
X 10~^ Z0 are the threshold for perturbation growth in 
these densities. However, we note that primordial gas 
in relic H II regions would be susceptible to ionization 
front instabilities at considerably lower metallicities be- 
cause they are compensated by larger free electron abun- 
dances. Our results indicate that this type of flow insta- 
bility would operate even in marginally enriched prot- 
galaxies and star-forming minihalos. The morphology of 
the instabilities at later times suggest that they might be 
efficient at driving turbulence in the shock, which could 
have important implications for the time scales of sub- 
seque nt star formation in the region (IWhalen fc NormanI 
l2007t l. Finally, we note that instabilities readily erupt 
at solar metallicities even with the much lower initial 
electron fr action Ye of 0-0001 than th e Xe = 0.01 value 
adopted in IWhalen fc NormanI (|2007f ). While in general 
it is difficult to know what free electron abundances are 
present in a galactic molecular cloud, a reasonable as- 
sumption is that all the carbon in the gas is singly ion- 
ized by the UV background. This yields a lower bound 




Fig. 9. — Dynamical instability in an r ^ density gradient in a variety of metallicities at 125.5 kyr. Panel (a): zero metallicity. Panels 
(b), (c), and (d) are for metallicities of 1 X 10"* Zq, 1 X 10"'^ Zq, and 1 X 10~^ Zq, respectively, while panels (e) and (f) are for 0.1 Zq 
and 1.0 Zq, respectively. 
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for electron fractions in star forming regions in the galaxy 
today that is on the order of those applied in this grid of 
metallicities, confirming that I-front instabilities will be 
a general feature of those regions. 

7. DISCUSSION AND CONCLUSIONS 

Instabilities readily form in primordial ionization 
fronts, even in the absence of metals. However, they 
emerge only when calculations couple multifrequency ra- 
diative transfer to gas phase primordial chemistry, which 
is crucial to capturing the formation of H2 that triggers 
dynamical instabilities in self-shielded shells or the high 
post-front gas temperatures that perpetuate shadow in- 
stabilities. We find that thin-shell modes cannot develop 
in atomic H and He but that Rayleigh- Taylor instabili- 
ties are not ruled out. Nevertheless, they do not appear 
on the time scales of UV breakout from the density pro- 
files in our models, which roughly approximate those in 
cosmological halos. 

Dynamical instabilities in D-type fronts were probably 
common in centrally-concentrated structures capable of 
molecular hydrogen self-shielding, like cosmological mini- 
halos at subparsec radii and molecular cores in primi- 
tive galaxies. However, shadow instabilities appearing 
earlier in R-type fronts and needing no H2 cooling may 
have pre-empted thin-shell instabilities in more evolved I- 
fronts. Inhomogeneities endemic to all primordial struc- 
tures made this instability ubiquitous, forming shocks 
and clumps that regulated the escape of ionizing UV pho- 
tons and metals from the first luminous objects. 

One question that remains to be answered is how insta- 
bilities forming at subparsec radii deep in cosmological 
halos or early galaxies governed the final morphologies of 
H n regions on kiloparsec scales. If clumping due to un- 
stable modes was long-lived, UV breakout along adjacent 
lines of sight may have been significantly delayed, causing 
the formation of extended fingers of ionized gas into the 
early IGM that resulted in serious departures from the 
classical " butterfly" appear ance of many cosmological H 
II regions (jAbel et al.ll2006( l. Numerical models bridging 
the spatial scales between the young D-type front and 
final H II region are now under development. 

Ionization front instabilities in the accretion envelopes 
of the first stars may have allowed them to grow to much 
larger masses than predicted by analytical models of in- 
fall reversal ([Tan & McKcc 2004). Fin gers of ionized gas 
elongating outward through the accretion shock might 
have channelled ionizing radiation out of the cloud that 
would otherwise have staunched infall. Accretion could 
then continue downward in columns of neutral flow, al- 
lowing the star to bypass the final mass limits set by the 
cutoff models. 

The strength and longevity of I-front instabilities de- 
termined their impact on subsequent star formation. 



Prominent extrusions from the front out into the halo 
or host galaxy might have created clumps of primordial 
gas prone to gravitational collapse by HD cooling in the 
relic H II region, forming a new generation of stars with- 
out recourse to chemical enrichment. This scenario mer- 
its further study with high resolution methods, as it is 
uncertain if clumps could have collapsed outside the po- 
tential well of the dark matter halo. On the other hand, 
more saturated modes due to less efficient cooling and 
fragmentation of the shell may instead have provided 
turbulent support against condensation in larger cloud 
complexes within protogalaxies, suppressing rather than 
promoting star formation. 

Clumping by instabilities may have determined the ex- 
tent of the first metal bubbles (| Johnson et al.|[2007f) in 
addition to modulating UV escape fractions and the fi- 
nal radii of H II regions. Mixing of metals with clumps 
may also have accelerated second generation star forma- 
tion. Efficiently cooled clumps would have collapsed to 
form stars on mass scales very different from their pre- 
decessors. Accurate models of primordial H II regions 
that include ionization front instabilities will establish 
the proper initial conditions for metal mixing in the next 
generation of early chemical enrichment studies. 

Although dynamical instabilities due to heavier ele- 
ments appeared at metallicities of 0.001 - 0.01 at ion- 
ized fractions of 0.01, they would have appeared in 
abundances an order of magnitude lower in the higher 
free electron fractions of fossil H II regions. This is 
comparable to the metallicity associated with the tran- 
sition from Pop HI stars to less massive populations 
(jMackev et al.ll2"003f ): hence, instability formation, and 
its consequences, smoothly followed rollovers in initial 
mass functions (IMFs) in the early universe. 

We neglect direct transport of recombination photons, 
which may limit the growth of unstable modes. Repro- 
cessed radiation emanating from within ionized exten- 
sions of gas may erode the sides of the fingers, photoe- 
vaporating overdensities in the troughs. This may not 
be an important effect in shadow instabilities since the 
interior of the jet is relatively rarefied, with low recom- 
bination emissivities. 
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